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Abstract. This paper aims to illustrate the applications of resonant Hamiltonian normal forms to some 
problems of galactic dynamics. We detail the construction of the 1:1 resonant normal form corresponding 
to a wide class of potentials with self-similar elliptical equi-potentials and apply it to investigate relevant 
features of the orbit structure of the system. We show how to compute the bifurcation of the main periodic 
orbits in the symmetry planes of a triaxial ellipsoid and in the meridional plane of an axi-symmetric 
spheroid and briefly discuss how to refine these results with higher-order approaches. 
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1 Description 



Nonlinear dynamics has often profited, since the pioneering works of Jeans, Lindblad, etc., of the progress in galactic 
dynamics. The interplay between the two disciplines has gradually established a common body of methods and results, 
as much as in the case of the interplay with celestial mechanics and is well documented in the literature [SlEltT] . 
Nowadays, advanced techniques of mathematical physics provide rigorous analysis of the main phenomena occurring 
in non-integrable dynamical systems potentially useful in galactic dynamics. However, the language of these modern 
methods often prevent a direct use of those results to applied fields like astrophysics. 

Aim of this work is to show how to apply modern perturbativc methods in analytical mechanics to provide a simple 
and comprehensive description of a common phenomenon in the dynamics of stellar systems: the onset of instability 
of a given motion due to a nonlinear resonance and the bifurcation of a new solution. We limit our analysis to the 
case in which the two involved frequencies, that of the perturbed solution and that of the 'perturbation', happen to 
become equal as functions in a suitable parameter space. This is called a '1:1 resonance' and plays a distinguished role 
in several concrete cases. The approach we use is that of mimicking the dynamics of the original physical system, that 
cannot be exactly solved, with an approximating system which is exactly solvable. The original system is assumed to 
be a generic non-integrable Hamiltonian system: the mimicking system is an integrable resonant Hamiltonian 'normal 
form'. We give explicit formulae to construct the normal form of a fairly generic class of systems and show how to use 
it to gather a general view of the phase space structure around the resonance. We prove that, in addition to a complete 
qualitative understanding, the theory also provides quite accurate quantitative predictions, even at the simplest level 
of the procedure. We finally discuss how to refine the predictions by going to higher-orders and how to generalize the 
method in a wider context. 

The plane of the paper is as follows: in section 2 we introduce the procedure to construct the approximating 
integrable system by recalling the method of the Lie transform |111I12| : in sections 3 and 4 we apply this approach to 
investigate some aspects of the dynamics in a symmetry plane of a triaxial ellipsoid pQ and in the meridional plane of 
an axisymmetric spheroid obtaining first-order estimates of the bifurcation thresholds of the 1:1 periodic orbits; 
in section 5 we discuss further developments and hints for other applications. 
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2 Normal forms 

The general setting in which we cast the dynamics of the physical systems we are interested in is that of a natural 
Hamiltonian in two degrees of freedom in which the potential is assumed to be expanded as a truncated series (in the 
coordinates q = (x,y)) of the type 

N k 

^feaj^^VilW^" 3 '. (!) 

k=0 j=0 

where the truncation order N and the parameters a are determined by the problem under study. We observe that 
for our 'galactic' applications N — 4 is already enough to describe the most relevant examples and moreover we will 
usually assume reflection symmetry in one or both the coordinate axes. The Hamiltonian 

H=\\p\ 2 + V^\q- a ), (2) 

is in the form of a power series and is therefore naturally apt to be treated in a perturbative way as a non-linear oscillator 
system. We construct a 'normal form' [5], namely a new Hamiltonian series which is an integrable approximation of 
the original one, suitably devised to catch its most relevant orbital features. 
The normal form is 'non-resonant' when the two harmonic frequencies 



u>i = ^2(7(2,0), w 2 = y 2(7(0,2)) (3) 

are generically non-commensurable: in this case we get explicit formulas for actions and frequencies of the box orbits 
parented by the x- and y-axis periodic orbits. A 'resonant' normal form is instead assembled by assuming from the 
start an integer value for the ratio of the harmonic frequencies and including in the new Hamiltonian terms depending 
on the corresponding resonant combination of the angles. This possibility might be thought to be exceptional: it is 
instead almost the rule because, even if the unperturbed system is non-resonant with a certain real value 

P = uj 1 /uj 2 (4) 

of the frequency ratio, the non-linear coupling between the degrees of freedom induced by the perturbation determines 
a 'passage through resonance' with a commensurability ratio, say mi/m^ with mi,m 2 € N, corresponding to the 
local ratio of oscillations in the two degrees of freedom. This in turn is responsible of the birth of new orbit families 
bifurcating from the normal modes or from lower-order resonances. 

We assume that our system is such that the ratio Q is not far from a rational value and then approximate it by 
introducing a small 'detuning' S so that 

p = mi/m 2 + S (5) 

and proceed like if the unperturbed harmonic part would be in exact m\\m2 resonance by putting the remaining part 
inside the 'perturbation'. We speak of a detuned (rax'-rri'i) resonance, implying with this a trick to deceive the formally 
correct but ineffective approach based on non-resonant generic frequency ratios. 

Resonant normal forms for the Hamiltonian system corresponding to ^ are constructed with standard methods: 
here we are going to apply the method based on the Lie transform [5] . The usual approach in Hamiltonian perturbation 
theory is to find a generating function in such a way to construct a 'simpler' Hamiltonian: the first successful algorithm 
devoted to a resonant system was precisely introduced to treat the problem of the 'third integral' in galactic dynamics 
|14) . The method based on the Lie transform offers several conceptual and technical advantages [MT2"] with respect to 
the original approach. The idea behind the method is that of seeing the canonical transformation as a 'flow' along the 
Hamiltonian vector field associated to the generating function. In spite of this apparently exotic statement, it turns 
out that this approach is that best suited in implementing recursive algorithm working with series of functions. In the 
following we provide the clue on how to set up the method remarking that for our applications only the first steps of 
the procedure will be necessary. 

Let us consider a phase-space function expanded as a power series in the canonical variables 

G = G + Gi + G 2 + ... (6) 

We assume that all functions are expanded in polynomial series and we adopt the convention that the 'zero order' 
terms of the series are quadratic homogeneous polynomials and terms of order n are polynomials of degree n + 2. 
To the series ^ is naturally associated the linear differential operator 
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whose action on a generic function F is given by the Poisson bracket: 

C G F = {F,G}. (8) 

After a scaling transformation 

Px — > y/uiPx, P v — > V^iPy, x — ► x/y/ui, V — > v/V^2, (9) 



the original Hamiltonian system ^ undergoes a canonical transformation to new variables (P, Q), such that the new 
Hamiltonian is 

K(P,Q)=c CG H(p,q), (10) 
where both functions are assumed to be in the form of power series 

N N 

H=J2 H n, K=J2 K n- (11) 

n=0 n=0 

To construct K starting from H is a recursive procedure exploiting an algorithm based on the Lie transform [5"lll2| . 
To understand how it works, let us consider the first step of the procedure, namely let us perform a transformation 
given by a function G\ considered as the first term in the generating function (we can start with the cubic term 
because Go gives only trivial linear transformations). The general relation (10 1 takes the form 

K + K 1 + ... = (1 + C Gl + ...)(H +H t + ...). (12) 

By equating polynomials of the same degree, we get the system: 

K = H , (13) 
K 1 = H 1 -£ Gl H , (14) 
K 2 = H 2 - L Gl H x - kC 2 Gl H (15) 

and equations involving terms of higher degrees. The first equality simply states that the zero-order new Hamiltonian, 

K = * (wi(Pf +X 2 )+u; 2 (P£ +Y 2 )), (16) 

coincides with the zero-order old (unperturbed) one. The second equation has to be solved to find the first order term 
K\. At this point we are faced with a difficulty: we have one differential equation involving two unknown functions, 
K\ and G\. To proceed we have to make some decision about the structure the new Hamiltonian must have, that is 
we have to chose a normal form for it. We therefore select the new Hamiltonian in such a way that it admit a new 
integral of motion, that is we take a certain function, say J, and impose that 

{K,I}=0. (17) 

The usual choice (but not necessarily the only possible) is that of taking 

I = H = K a (18) 



so that the function ( 16 ) plays the double role of determining the specific form of the transformation and assuming 
the status of the secona integral of motion. With this choice, the fundamental equation of the chain, that we can also 
write in the form 

K 1 = H l -L Gl H = H l + L Ho G l , (19) 

is solved with a trick that we illustrate in the following. The action of the operator Ch on any polynomial function 
which commutes with H$ (and therefore that has Hq as an integral of motion) is to kill it, whereas its action on 
any other polynomial gives a uniquely defined non-vanishing polynomial. We can therefore split the polynomial Hi 
appearing in ( 19 ) according to the rule 

Hx = Hf + Hf (20) 
where H^ is the part which stays in the kernel of Ch , that is is such that 

C Ho Hf = 0, (21) 
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and Hf~ is the part which stays in the range of Ch , that is is such that 

Ch H^ = Ri, (22) 

where Ri is a non-vanishing cubic polynomial. Since our new Hamiltonian, with the choice made, is in normal form if 
and only if it stays in the kernel of £h , we can then solve eq.(19) by applying the simple prescription: 

K x =Hf, G l = C H \Hf. (23) 

We have then constructed the normal form Kq + K\ at order 1 and computed the first term of the generating function 
G\. we can therefore use it in the subsequent equations of the system to compute the terms K n ,n > 1 (which are 
still not in normal form) and go one step further by expanding eq.(12| at order two and repeating the procedure to 
compute G 2 and the normal form at order 2 and so forth. 

Formally, a more direct way of applying this method is by using 'action-angle'-Zifce variables J, 6 defined through 
the transformation 

X = v / 2^lcos0i, F = V^2Cos6>2 ! (24) 
P x = v/2Jism(9i, P Y = yj%hsm.6 2 . (25) 



In this way we have 
so that 



Hq = u>i J\ + w 2 J 2 , (26) 

Q Q 

£^=^—+^2 — . (27) 

Moreover, a generic polynomial series turns out to be a Fourier series in the angles with coefficients depending on the 
actions. It is clear that in the non-resonant case, the kernel of Ch is composed only of functions of the actions and 
this must also be the case for the normal form: in this case we speak of a Birkhoff (or non-resonant) normal form. But 
if the frequencies are (almost) in the ratio 77*4 /m 2 , the typical structure of the resonant normal form (truncated when 
the first resonant term appears) is [7] 

m 1 +m 2 



K = m 1 Ji+m 2 {J 2 + SJi)+ J2 l ?{k) (Ji,J 2 ) + 



k=2 



*mim2 



J™ 2 J™ 1 cos[2(m 2 6<i - nMa)], (28) 



where V {k) are homogeneous polynomials of degree k whose coefficients may depend on 5 and the constant a mim2 (a) 
is the only marker of the resonance. It is easy to check that this is the most general form of a phase-space function of 
degree m\ + m 2 in the actions which stays in the kernel of Cr as given by (27) if the frequencies are such that, in 
agreement with assumption ([5]), 

mi/m2 = u>i/lo 2 + 8. (29) 

In these variables, the second integral is 

£ = m\J\ + m 2 J 2 (30) 
and the angles appear only in the resonant combination 

if; = m 2 0i - mi9 2 . (31) 

For a given resonance, these two statements remain true for arbitrary N > mi Introducing the variable conjugate 

to ip, 

1Z = m 2 Ji — mi J 2 , (32) 

the new Hamiltonian can be expressed in the reduced form K(TZ,ip;£^a), that is a family of 1-dof systems in the 
phase-plane 1Z, ip, parametrized by £ (and a). 

We recall that, as in ([5]), the resonant ratio is not exact: as a result we have that the quadratic part can be written 
in the form 

#o = ^2 (—Ji + J2) + U2SJ1. (33) 

We can then rescale the Hamiltonian with the constant factor m 2 /w 2 and normalize with respect to the exactly 
resonant quadratic part (30): the small 'detuning' term m 2 (5Ji, even if still quadratic in the coordinates, is treated as a 
higher-order term and inserted into the perturbation H\ + H2 + ... Since in practice 6 is related to some morphological 
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parameters of the potential, this trick allows us to considerable extend the range of the analysis. In the applications 
we aim at illustrating the fundamental case in which the ratio between the integers is simply one. Therefore, we will 
explicitly compute the 1:1 resonant normal form that, when truncated to the first non zero resonant term, is given by 

K = J x + J 2 + 5J X + aJ\ + bJl + ,hJ 2 (2c + dcos[2(0 a - 6 2 )}) , (34) 

where the coefficients a, b, c and d = an are the result of the normalization process and are completely defined in 
terms of the parameters of the original system. 

In the applications below we arc interested in the global structure of phase-space, but also the explicit solution 
of the equations of motion is of great relevance. For a non-resonant normal form, the problem is easily solved: the 
coefficient a mim2 vanishes and K no longer has terms containing angles. Therefore, the J are 'true' conserved actions 
and the solutions are 

X(t) = v/2Ji"cosKir, y(r) = v/2J2C0s(k 2 t + 6» ), (35) 

where 

k = VjK (36) 

is the frequency vector and 8q is a suitable phase shift. In the resonant case, in general it is not possible to write the 
solutions in closed form. It is true that the dynamics described by the 1-dof Hamiltonian K(TZ, ip) are always integrable, 
but the solutions cannot be written in terms of elementary functions. However, solutions can still be written down in 
the case of the main periodic orbits, for which J, 9 are true action-angle variables. There are two types of periodic 
orbits that can be easily identified: 

1. The normal modes, for which one of the J vanishes. 

2. The periodic orbits in general position, namely those solutions characterized by a fixed relation between the two 
angles, m 2 0i — mi6 2 = 8 . 



In both cases, it is straightforward to check that the solutions retain a form analogous to Eq. (|35j) with known expressions 
of the actions and frequencies in terms of £ and the parameters a such that k\/k 2 = mi/m 2 . By using the generating 
function Eq.([6]), the solutions in terms of standard 'physical' coordinates can be recovered (taking into account possible 
scaling factors) by inverting the canonical transformation. The transformation back to the physical coordinates gives 
a series of the form 

Q(t) = 9i + 92 + 93 + ■■■ (37) 
and can be computed by exploiting the properties of the Lie series as above: 

91 = Q, (38) 

92 ={Gx,Q}, (39) 

93 = {G 2 ,Q} + i{Gi,{G 1 ,Q}} (40) 



and so forth. From a knowledge of the normalized solutions Eq.(35l, we can therefore construct [3] power series 



approximate solutions of the equations of motion of the original system 

^ = -V.vW(x,y;a). (41) 

As a rule, normal modes exist on each 'energy' surface K = E = (m 2 /uj 2 )E. Periodic orbits in general position 
exist instead only beyond a certain threshold and we speak of a bifurcation ensuing from a detuned resonance. The 
bifurcation is usually described by a series expansion of the form 

E = J2*k\6\ k , (42) 

k 

where the are coefficients depending on the order N and the parameters a.. The order of truncation of the series 



is related to that of the normal form. Eq. ( 42 ) implies that at exact resonance (vanishing detuning) the bifurcation 
is intrinsic in the system and that, going away from the exact ratio of unperturbed frequencies, gradually increases 
the threshold value for the bifurcation. We will see that already a linear relation given by the first order truncation 
provides a reliable estimate of the threshold values if due care is made of the consistency of series expansions. 
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3 The stability of axial orbits in systems with elliptical equipotentials 

We show how to exploit the above theory in some cases of interest in galactic dynamics. We start with the problem of 
the stability of axial orbits in triaxial potentials and of the possible existence of additional stable families of periodic 
orbits. The relevance of such issues in galactic dynamics is for problems like the construction of self-consistent equilibria, 
the computation of isophotal shapes and velocity ellipsoids, etc. 

We consider a fairly general class of potentials with self-similar elliptical equipotential and unit 'core' radius of the 
form 2113 

a I ?J » r (43) 

|log(l+x 2 + ^), a = 0. 

The ellipticity of the equipotentials is determined by the parameter q: we have an 'oblate' figure when q < 1 and 
a 'prolate' figure when q > 1. The range of q is constrained by the condition that the density contours of the mass 
distribution generating the potential resembles real galaxies. If the slope parameter is restricted to the range 

- 1< a < 2, (44) 

Evans [9 shows that the allowed range for q is an interval around 1 whose extent depends on a and shrinks to zero if 
a tends to —1. For example, in the logarithmic case (a = 0) the range to get bulges or ellipticals is 

0.7 < q < 1.1 (45) 

whereas, for a tending to +1 the range tends to the interval 

< q < 1.25. (46) 

The family of potentials ( 43 ) can be expanded in a series of the form ([I]) . In view of the reflection symmetry with 
respect to both axes, the non-vanishing coefficients of the expansion up to order N — 4 are the following: 

(47) 
(48) 
(49) 
(50) 
(51) 



C(2,0) — 


1, ,2 _ 1 
2 W 1 - 2' 


C(0,2) = 


1, ,2 _ 1 
2 W 2 - 2 9 2 


C(4,0) = 


a-2 
8 ' 


C(2,2) = 


a-2 
4q 2 » 


C(0,4) = 


a-2 
8q 4 ■ 



The normal form truncated to the first non zero resonant term is given by expression (34) where, as can be verified 
by using the procedure above, 6 = q — 1 and the coefficients are given by 

~- 3 ° {m ~\ q (*-2), (52) 



4 C"(2,o)\/2C(o,2) 16 

b= Jr /or = Tfi^ (Q ~ 2) ' (53) 

4 ^(0,2) V 2C/ (0,2) Lb 1 

, 1 C(2.2) 1, „s ,_.s 

c = d =Jr /or — = s (a ~ 2) ' (54) 

4 C (0,2) V^^.O) 8 



The model system given by ( 43 ) represents motion in the symmetry planes of a triaxial galaxy. In each of those 
planes, the symmetry axes directly give periodic orbits: at low energy, since the dynamics are slightly different from 
those of a harmonic oscillator (a = 2), we may expect them to be stable oscillations. What happens when energy 
increases? Nonlinear dynamics give asynchronous motions, with frequencies depending on amplitudes. Instability can 
be triggered by low-order resonance and we can see how the normal form ( 34 1 is ideally suited to investigate this 
phenomenon. 

By introducing the canonical variables adapted to the resonance by the resonant combinations 

■R = ,h-J 2 , ^ = 2(0i-0 2 ), (55) 
the new Hamiltonian can be expressed in the reduced form 

K = £ + -5(8 + TZ) + A(£ 2 + TZ 2 ) + B£TZ + C(£ 2 - K 2 ){2 + cos V), (56) 
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with 

A=\(a + b), B= l -{a-b), C =| (57) 

and 

£ = Ji + -h (58) 
is the second integral as in (30 1. Considering the dynamics at a fixed value off, we have that K defines a one-degree 



of freedom (ip, TZ) system with the following equations of motion 

ip = ~5 + B£ + 2(A-C(2 + cosip))lZ, (59) 

1l = C{£ 2 -K 2 )smip. (60) 

The fixed points of this system give the periodic orbits of the original system. The pair of fixed points with TZ = ±£ 
correspond to the normal modes: the orbit TZ = £ (namely J2 = 0) is the periodic orbit along the x-axis (long axial 
orbit if q < 1) whereas the orbit TZ — —£ (Ji = 0), is the periodic orbit along the y-axis (short axial orbit if q < 1). 

However, we can see that additional periodic orbits may appear. These periodic orbits 'in general position' exist 
only above a given threshold when the axial orbits change their stability. This phenomenon can be seen as a bifurcation 
of the new family from the normal mode when it enters in 1:1 resonance with a normal perturbation. The phase between 
the two oscillations also plays a role: these additional periodic orbits are respectively given by the conditions ip = 
(inclined orbits) and ip — ±tt (loop orbits). These conditions are solutions of TZ = (when TZ 7^ ±£) and determine 
the corresponding solutions of ip = 0: 

5 + 2B£ 5+(a-b)£ 
Ui =4(3C^A) = ic-a-b (61) 

and 

5 + 2B£ S + (a-b)£ 

^ = ± ^ n ?=A(C^A) = c-a-b • 

In view of (55 1 and (58), the constraints < J%, J2 < £ applied to these solutions translate into the conditions of 
existence 

S S 

£>£a = o 7T or £> £ i2 = — — — (63) 

3c — 2a 2b — 3c 

and 

£>£ a = or £>£&=- . (64) 

~ c-2a ~ 2b -c K ' 

Index I denotes the critical values at which the loop orbits bifurcate (ther e ar e two of them, one clockwise, the other 



counter-clockwise) respectively from the rc-axial normal mode if the first of ( 64 ) is satisfied and from the y-axial normal 



mode if the second one is satisfied. Index i denotes the critical values at which the inclined orbits bifurcate, respectively 



from the x-axial normal mode if the first of ( 63 1 is satisfied and from the y-axial normal mode if the second one is 
satisfied. At the same bifurcation values, the normal mode suffers a change of stability, passing from stable to unstable 
when the new orbit is born. Also a second transition is possible (unstable to stable) if a second bifurcation ensues. 

In fig [I] we can see the dynamics of the 1:1 doubly-symmetric resonant Hamiltonian displayed with the two coor- 
dinate sets exploited above: the left panels contain the level lines of the reduced Hamiltonian in the (ip, TZ) cylinder 
with the lines ip = and ip = 2ir identified; the right panels are the (X,Px) 'surfaces of section' computed with the 
condition Y = 0, Py > 0. We have chosen a model with an oblate (q = 0.7) 'logarithmic' potential (a = 0). The upper 
panels are computed for £ = 0.3, that is below £ti(— 0.365): on the left we see the continuous curves associated to the 
'rotational' invariant tori around the two stable normal modes; on the right we see a section of these tori corresponding 
to 'box' orbits around the minor-axis. The lower panels are computed for £ = 0.7, that is above £ii'- on the left we see 
that at ip — 7r a stable fixed point has appeared surrounded by 'librational' invariant tori; on the right we see that the 
minor-axis is now unstable and there are two stable fixed points corresponding to the two periodic loops orbiting in 
the two senses of rotation. Each of them is surrounded by non-periodic loops. There is also a separatrix beyond which 
there are the remaining boxes parented by the mayor-axis which is still stable. 

Fig{T] nicely provides qualitative informations. We stress that these results are also quantitatively useful. In view of 
the rescaling and of the expansion of the energy as a truncated series in the parameter £, we have that E — u)2£ = £/q 
is a first order estimate of the 'true' energy of the orbital motion. We can use the above critical values to establish 



the instability threshold for the model problem given by potentials (43). For the parameter ranges used in galactic 



dynamics, it happens that the typical bifurcation is that of the loop orbits at the critical energy 



(65) 
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Fig. 1. Dynamics of the 1:1 doubly-symmetric resonant Hamiltonian, see text for explanation. 



In the range 



0.7 < q < 1.3, 



(66) 



which can be considered as 'realistic' for elliptical galaxies, the thresholds (65) give estimates correct within a 10% if 
compared to numerical computations P0I2- 



4 Resonances and bifurcations in axisym metric potentials 



Here we discuss how a slight generalization of the normal form studied above is able to describe some very important 
features of the dynamics of axisymmetric galaxies. A spheroidal bulge or a disk/halo system can be modeled with 
a potential with an exact axial symmetry and the conservation of the angular momentum L with respect to the 
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symmetry axis still allows us to stay with only two degrees of freedom. The dynamics is then given by the family of 
potentials 

L 2 

V cS (R, z; L, q, a) = ^—^ + V(R, z; q, a), (67) 



where potentials of the same form as in (43) can be used by exploiting cylindrical coordinates R, (j), z. These effective 
potentials have a unique absolute minimum in 

R = R c (a) = L^, z = 0. (68) 

This is a stable equilibrium on any meridional plane <f> = const corresponding to a circular orbit of radius R c (a) of 
the full three-dimensional problem. Now we have four parameters, L, q, a and the energy to take into account. We can 
get rid off of the energy by putting to zero the core radius so that the dynamics are scale-free and we may fix energy 
at some convenient value: 

(! + i)e-*fe, a^O, 
0, a = 0. 



E — E a - 

This implies that the radius of the circular orbit at this energy is 

R c (a) = e~5T^, -l < a < 2 
and we can investigate the dynamics at E = E a by varying L in the range 

< L < L max = i?, 



(69) 



1 



(70) 



(71) 



without any loss of generality [T51I4] . 

In order to implement the perturbation method we expand the effective potential around the minimum ( 68 ) . We 
introduce rescaled coordinates according to 



R — R c 

Rc. 



Rc 



(72) 



with origin in the equilibrium point (68). The potential (67) is expanded as a truncated series (in the coordinates x, y) 
of the form ([!]) , where the truncation order N is determined by the resonance under study: here we limit ourselves to 
the 1:1 case. From (68) and the rescaling (72), the constant term of the expansion is 



L, a = 



C(o,o)(L, a) = 
and the other coefficients have the form 

C(j,k-j)(L,a,q) = L^sc {j>k _j){a,q). 



In order to simplify formulas, we introduce the new parameter 

2a 



2 + a' 



K P < 2, 



(73) 



(74) 



(75) 



with the same range of a in view of ( 44 ) . 

The orbit structure of the original family of potentials ( 67 ) at the energy level fixed by ( 69 ) will be approximated 
by the orbit structure of the rescaled Hamiltonian 



where 



H = +P 2 y ) + Ves{x,y;q,a), 



V cS {x,y;q,a) = L^Veg (x, y; L, q, a). 



The dynamics given by Hamiltonian ( 76 ) take place in the rescaled time 

r = t/L^ 1 



(76) 
(77) 

(78) 
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at the new 'energy' 

E F = (E a - C m (L, a)) = 1 (l - (L/L max f) . (79) 



According to ( 71 ) the singular value L = is excluded from the analysis implying that the fictitious energy ([79| is 



always finite and that the expansion around the equilibrium point (68) make sense. The equilibrium value Ep = 
corresponds to the circular orbit with L = L max and lower values of L give increasing values of Ep. 
The non-vanishing coefficients of the expansion of V e g up to order TV = 4 are the following: 

c (2 ,o) = |w? = (80) 

C(0,2) = \"\ = , (81) 

10+3a-q 2 {qo\ 

C(3,0) = g ' ( 82 ) 

C (l,2) = ( 83 ) 

n, , — 54+lla-6a 2 +a 3 foA\ 

C(4,0) — 24 ' \ m > 

C (2,2) = 4^2 , (85) 

C(o,4) = ~ispr- ( 86 ) 



The first two of them provide the frequencies of the epicyclic motions. Recalling (|3| and the time rescaling in ( 78 ) 
the radial and vertical harmonic frequencies are respectively 



L/3+1 L/3+1 

and 

uj 2 1 



V2+^ (g7) 



L 0+l qL fl+l ■ 

Therefore, to describe the bifurcations ensuing from the 1:1 resonances, the normal form is computed by using the 
detuning parameter defined through 

6 = - - 1 = q\/2 + a - 1. (89) 



v 



The presence of coefficients of cubic terms in the expansion accounts for the breaking of the reflection symmetry 
with respect to the y axis: we still have the reflection symmetry with respect to the equatorial plane corresponding to 



the x-axis in the coordinates (72 ). The X-axis ( J\) normal mode corresponds, in the full 3-dimcnsional problem, to the 
equatorial non-periodic disk orbits. Breaking the symmetry with respect to the 'vertical axis' means that the y-normal 
mode corresponds to another orbit that in the 3-dimensional space is confined to a folded surface: the so called thin 
tube. This is an example in which we see that normalizing transformation are, in a sense, 'rectifying' transformations. 
On this respect, it is also useful to recall something about the taxonomy of periodic orbits. Again we denote the 
bifurcating family as the inclined orbit in view of its natural interpretation as the in phase (tp = 0) 1:1 resonance of 
the two oscillations [TUTU]. The anti-phase resonant loops never appear in these systems, at least as a stable family (see 
below). The inclined periodic orbits parent two families of inclined boxes that may arrive quite far from the equatorial 
plane both above and below the disk: this phenomenon is called levitation [17] . We recall that our inclined orbits have 
also been referred to as reflected banana by 18J and [SJ and simply as banana by |15j : we prefer to leave this term 
as the standard denomination [19] for the 2:1 resonance. They exist also nicknames for high-resonant periodic orbits 
(fish, pretzel, etc. [T51IT5] ). 



The normal form truncated to the first non zero resonant term is given by the general expression ( 34 ) where the 
coefficients are 

3 c (4 ,o) 15 45,0) _ ?(-5g(5 - a) 2 (2 + a) + 3\/2T^(27 - 8a + a 2 )) 

, (911) 



2 lo 2 lo 2 4 ufu% 48a/2 ■ 

3c (0 ,4) 5 cfi, 2 ) 1 , /9 , 5(2 -a) 



a 



b="^l-~ ^ = — (a -2) , (91) 

C( 2 ,2) 3c (li2 )C ( 3 j0) c 2 {l 2) (io - 7 a + a 2 )(l - 3q^2 + a) 

! ( yz ) 



2wiw| 2w^w| 3wiw| 24^2 + a 

d = C( 2 .2) C(i, 2 )C(3,o) _ C(i, 2 ) = (a - 2) (3 - 5qV2 + a - a(3 - qy/2 + a) 
_ 2u} X ujj 2wfw| ~ 24V2 + a 
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Comparing with eq.(54) we see that the effect of the symmetry breaking is that now d ^ c with more complex 



bifurcation conditions. 

By lowering the angular momentum (namely by increasing the fictitious energy Ep) both normal modes may lose 
their stability through a 1:1 resonance. The equations of motion for the variables adapted to the resonance now are 



ip = K n = ^5 + B£ + -(4A - (2c + dcostp))n, 
K = -K 4 , = j (S 2 - K 2 ) sinV>. 



(94) 
(95) 



The analysis of the fixed points of this system proceed as in the previous section. The pair of fixed points with 1Z = ±£ 
correspond to the normal modes. The periodic orbits in 'general position' (inclined and loops) are respectively given 
by the solutions: 

6 + 2B£ 5+ (a- b)£ 



and 



ip = 



1p = ±7T 



Hi 



2c + d- 4A 2c - a-b + d 
S + 2B£ 5 + (a- b)£ 



2c-d-4A 2c-a-b-d 
The constraints < J\ , J2 < £ applied to these solutions translate into the conditions of existence 

S „_ „ S 

~~ ~d 



£ - £il = ^7 TT^ or £ - £i2 = ^7T \ 

2(c — a) + d 2{b — c) 



and 



£>£ 



11 



or £ > £12 = 



2(6 



(96) 
(97) 

(98) 
(99) 



By using (90 93) the bifurcation equations ( 98|99 ) determine the critical values of £ in terms of the parameters q, a. 
A common situation, for models ranging from sensibly oblate to prolate and with a in the range ( 44 ) , is that in which 



the IE-axis becomes unstable and the inclined appears as a bifurcation from the disk orbit. The passage to instability 
of the y-axis is possible for oblate models such that 



Q < 



1 



V2- 



(100) 



and gives rise to a bifurcation from the thin tube. In this case, also the second bifurcation is usually possible: the 
disk orbit returns stable and the new family (loops) is unstable. In fig{2]we can see the dynamics of the 1:1 resonant 
Hamiltonian with one of the symmetries broken, displayed with the same settings as in fig[l] We ha ve ch osen a model 
with an oblate 'logarithmic' potential (a — 0) with a value of q = 0.65 in order to satisfy inequality ( 100) to show the 
occurrence of both bifurcations. The upper panels are computed for £ = 0.13, that is above £ii{= 0.10): on the left 
we see the island at ip = 0, 2n around the stable inclined orbits bifurcated from the thin tube; on the right we see a 
section of the tori around the stable inclined orbit, a separatrix departing from the thin tube (the unstable normal 
mode) and the tori beyond it associated to the box orbits parented by the stable disk. The lower panels are computed 
at £ = 0.15, that is above £n{= 0.135), the threshold of the second bifurcation: on the left we see that at ip = ir an 
unstable fixed point has appeared, testifying the instability of the loops ensuing from this second bifurcation; on the 
right we see that the the thin tube is again stable and there are two unstable fixed points corresponding to the two 
loops. The thin tube is now surrounded by tori associated to the thick tubes and there the other two families of tori 
around the inclined and around the disk. 

To make quantitative predictions, we want expressions for the corresponding critical angular momentum. The 
approach we have followed so far is altogether a perturbation approach truncated to the first non-trivial order: therefore 
it is natural to look also in this case for e xpan sions truncated to the first o rder in the detuning parameter. Taking 
into account the rescaling in K, conditions ( 98 1 and the explicit expressions ( 90 93 ) , the first order expansions of the 



critical values of the fictitious energy ( 79 ) are 



12(2+0) 



Ep = 



<5, S < 0, 



5(-2-q+q 2 
2-\-a — a z ' 



(101) 



The first solution corresponds to the bifurcation from the thin tube, the second one corresponds to the bifurcation 
from the disk. These are again examples of series of the form (42) truncated to the first order. Afterwards, using the 
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Fig. 2. Dynamics of the 1:1 simply-symmetric resonant Hamiltonian, see text for explanation. 



relation between Ep and L established by ( 79 1 , we get the following expressions for the critical values of the angular 
momentum below which inclined orbits exist: 



L, 



crit — i— x 



1 24a(gy^+a-l) ^ ^ < 1 



5(2+a-a 2 ) 

120(9^2+0-1) 
2+a-a 2 



, q > 



2+a ' 
1 



(102) 



/2+a ' 



It is also useful to write the limiting case of the logarithmic potential (a = 0): 



29 i 12 /n 1 

e io+ 5 V2 9j g < 



(103) 
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A comparison with the outcome of numerical determinations of the bifurcation threshold allows us to evaluate the 
accuracy of these analytical predictions. In [3] the critical value of the angular momentum for the bifurcation of the 
inclined orbits, computed with Eq.( 102[ ) for general a and with eq.(103) for a = 0, are compared with numerical data 



obtained either from published works 18,9,15 or by numerical computations made specifically for that paper. In this 
case, the bifurcation has been detected tracing the instability threshold of the normal mode by means of the Floquet 
method. The accuracy is particularly good when the model is close to the exact resonance. Overall, the discrepancy 
linearly grows with detuning, as can be expected in this first order approach. In the case of two strongly oblate models, 
the thin tube becomes unstable: in the rather extreme case with q = 0.4, a = 0.5, that is just at the margin of the 
physical range established in [9], the detuning is 6 — —0.37 and the relative error in the prediction is 18%. In all other 
cases, the disk becomes unstable and the quality of the prediction can be represented by the cases with q = 0.8, a = 0.1 
studied by Hunter and collaborators |15| . where the detuning is S — 0.16 and the relative error in the prediction is 8% 
and q = 0.85, a — —0.18 which was investigated by Evans [9], where the detuning is S = 0.14 and the relative error in 
the prediction is 7%. 

As said above, a second stability change may occur when the second resonant family appears [T5]. In this setting, 



this possibility occurs if the loops appear. By using the bifurcation equations (99 1 and the explicit expressions (90 - 93 ), 
the inequality can be satisfied only with values of the parameters corresponding to rather extreme oblate models. This 
implies a negative value of the detuning and a critical fictitious energy 



2 + a — or 



We get the following expression for the critical value of the angular momentum below which loops bifurcate from the 
thin tube: 



1 80^/2+^-1) _ 



Ve \ 2 + a 

This is again a 'pitchfork' bifurcation so that the thin tube regains its stability. The loops are unstable and, lowering 
the angular momentum below the critical value, remain unstable for every reasonable combinations of the parameters. 
At the same time, the inclined tend to occupy an even larger fraction of phase space. We are now very far from the 
harmonic behaviour, therefore it is reasonable to expect a certain discrepancy between the prediction and the actual 
value. Referring again to the models in [J], in the case with q = 0.4, a = 0.5 the 'true' critical value for the return to 



stability is 0.32 [TS] whereas (105) predicts 0.18. In the concluding remarks we give some hints on how to improve the 
quality of the predictions. 



5 Discussion 

Our aim has been to show how to make useful predictions on the dynamics with the 'tool-box' offered by the normal- 
form theory. The reader can also simply choose one of these tools and use it in some specific problem (if it can be 
fitted to the classes above) even without a deep involvement into the theory. However, we have also tried to lay out 
a sketch of how to proceed from the oustart. We remark that this approach is not the only possible one and that 
many others have been devised: after all, many of the phenomena described above are a form of nonlinear resonance 
in mechanics, a subject to which many efforts have been devoted [20] in the long history of this field. However, to 
testify the relevance of the normal-form method, it has also been applied as a tool to investigate parametric resonance 
in its more general mathematical formulation |21j . 

In order to improve the predictive power of the approach, the standard way to proceed is that of going to higher 
orders in the expansions. Following the procedure illustrated in section 2, an algorithm can be devised to construct a 
normal form up to arbitrary order. With this, one can compute the bifurcation thresholds as series in the detuning 
truncated at an order to be determined by the required precision. Very good predictions can be performed with series 
truncated already at degree 6-8 [T][2] m the phase-space variables for several kinds of resonance. Clearly, an immediate 
warning to do is that the procedure may result in cumbersome formulas, therefore some educated exploitations of 
automatic computations is mandatory. But there is also a more subtle issue with high-order expansions: the series 
we are working with are in general not converging. However, they are nonetheless asymptotic series. This means that 
there exist an 'optimal' order at which to stop anyway the expansions because beyond it nothing is gained in precision 
(for a discussion see ref.[2]). There is in general no way to know the optimal order in advance, however experience with 
known systems and comparison with numerical analysis suggests that, when/if the system enters a chaotic regime (e.g. 
when energy overcomes a certain stochasticity threshold |22j ) the optimal order apt to describe the dynamics near to 
the chaotic transition is just the order at which resonant terms firstly appear in the normal form. 

As far as the applications in galactic dynamics, there are several instances in which some progresses can most 
probably be done with resonant normal forms. We mention the study of rotating systems in which not much more 
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than the results obtained in the pioneering works of de Zeeuw and Merritt |10| is available from the analytical side 
and the investigation has recently be done only numerically. Another issue of great relevance is that of a central black 
hole (or some other kind of gravitational singularity like in a strong galactic cusp) in which the dynamics are that of a 
quasi-Keplerian perturbed field. Regularizing coordinate transformations allows us to turn the system into a perturbed 
oscillator and the normal form method can again be applied. Finally, it is clear that an important step forward is 
that of investigating fully tri-dimensional systems. In this case, in general, the normalizing transformation no longer 
provides an integrable normal form. The gain is usually limited to decrease by one the dimensionality of the problem. 
However, the investigation of existence and stability of periodic orbits proceeds in much the same way as above, even 
if at the price of a higher algebraic complexity. This is probably the most promising area of future investigations. 
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